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There has been a lot of research interest in modified gravity theories which utilise the Vainshtein mechanism to 
recover standard general relativity in regions with high matter density, such as the Dvali-Gabadadze-Porrati and 
Galileon models. The strong nonlinearity in the field equations of these theories implies that accurate theoretical 
predictions could only be made using high-resolution cosmological simulations. Previously, such simulations 
were usually done on regular meshes, which limits both their performance and the accuracy. In this paper, we 
report the development of a new algorithm and code, based on ECOSMOG, that uses adaptive mesh refinements 
to improve the efficiency and precision in simulating the models with Vainshtein mechanism. We have made 
various code tests against the numerical reliability, and found consistency with previous simulations. We also 
studied the velocity field in the self-accelerating branch of the DGP model. The code, parallelised using MPI, is 
suitable for large cosmological simulations of Galileon-type modified gravity theories. 



I. INTRODUCTION 

Modified gravity theories 0]] are proposed as alternative to 
the dark energy scenarios to explain the observed acceler- 
ating expansion of our Universe yHil], and have attracted a lot 
of research interest recently. Instead of invoking some sort of 
mysterious new energy component which drives the dynamics 
of the cosmos, these theories suggest that the Universe is filled 
with only normal plus dark matter (which is usually assumed 
to be cold), but the law of gravity can be different from that of 
standard general relativity (GR) on large scales, leading to the 
speed-up of the expansion rate. 

Since the law of gravity is considered as universal, a modifi- 
cation to GR on large scales almost necessarily implies corre- 
sponding changes in the behaviour on small scales. Any such 
changes from GR, however, is already highly constrained by 
numerous local tests of gravity H, and this has rendered many 
models incompatible with experimental bounds. Thus, any vi- 
able modified gravity theory should have some mechanism by 
which such modifications are suppressed and GR is restored 
in high-density regions such as the Solar system, where those 
experiments have been carried out and the resulted bounds ap- 
ply. Such mechanisms are commonly referred to as 'screening 
mechanisms' in the literature. Obviously, to make the model 
theoretically appealing, the screening mechanism needs to be 
an inherent (rather than add-on) property of it, which comes 
from the dynamics of the theory. The screening effect means 
that gravity behaves in different manners in different environ- 
ments, and such environmental dependence often boils down 
to a high degree of nonlinearity in the relevant field equations, 
which makes the study of such theories rather challenging and 
(numerically) demanding. 

In most modified gravity theories studied so far, the modi- 
fication to GR is realised by a dynamical scalar-type (spin-0) 



* Email address: baojiu.li@durham.ac.uk 
t Email address: gong-bo.zhao@port.ac.uk 
t Email address: kazuya.koyama@port.ac.uk 



field, which mediates the modified gravitational force. These 
theories could be roughly divided into two classes. In the first 
class of such theories, the screening is achieved by nonlinear 
couplings of the scalar field to matter and/or a nonlinear po- 
tential governing the self-interaction of the scalar field. If the 
coupling and potential are appropriately specified, the scalar 
field may acquire a very heavy mass in dense regions, making 
it hardly propagating or mediating any force between matter 
particles, or extremely weakly coupled to matter such that the 
resulted modification to GR very small. The chameleon mod- 
els EKH, with f(R) gravity model d (see also (THH) 
as a special example, belong to the former case, while the dila- 
ton 111611 and symmetron 01711 models belong to the latter case. 
In these models, f(R) gravity (mainly the model of J3) has 
been the most well-studied, and there have been many works 
which studied in detail its structure formation in the nonlin- 
ear regime, with the aid of TV-body simulations lfl8l - [3lll . Dur- 
ing this process, we have developed an efficient TV-body code, 
ECOSMOG Ull], based on the publicly available code RAMSES 
ll33ll . which is massively parallelised using MPI and therefore 
makes large simulations for f(R) gravity possible. Using the 
generic parameterisation for modified gravity theories of the 
first class JMIHl, ECOSMOG has been recently extended to 
simulate general chameleon, dilaton and symmetron cosmolo- 
gies sin. 

The other class of modified gravity theories involving scalar 
degrees of freedom, with the Dvali-Gabadadze-Porrati (DGP) 
brane-world model ll38ll as the most well-known example, re- 
alises the screening using nonlinear derivative self-couplings 
of a scalar field. Here, the nonlinearity makes the gravitational 
fields of individual particle interfere so that the deviation from 
GR per particle is weakened in a collection of particles com- 
pared with that for an isolated particle. This is known as the 
Vainshtein mechanism ll39ll in the literature, which was origi- 
nally introduced in the context of massive gravity to suppress 
the extra helicity modes of massive graviton to recover GR in 
the massless limit. The Vainshtein mechanism occurs not only 
in the recently proposed nonlinear massive gravity IHOWa 
and braneworld models but also in more general setups, such 
as the Galileon model ll43l - l48ll . The cosmological implications 
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of the Galileon model have been studied by many groups (see, 
e.g.. ll49l - l60ll ) in details. However, to date studies of large-scale 
structure formation in the nonlinear regime, even for the sim- 
plest case of DGP, have been very limited. The primary reason 
is that the field equations for these theories generally involve 
high products of second-order derivatives (see below), which 
make them difficult to solve numerically. There have been a 
few works on iV-body simulations for DGP ll6lU64ll . but these 
are so far limited to codes which, unlike ECOSMOG, have no 
mesh refinements. Furthermore, for the code to solve the DGP 
equation one may have to smooth the underlying density field 
to reduce the noise, which certainly limits the force resolution 
of the simulations. 

In this paper, we report the development and initial results 
of a new iV-body code, ECOSMOG- V, which is a sibling ver- 
sion of ECOSMOG but for simulations of theories involving the 
Vainshtein screening. We have made various tests of the code 
to be confident about its accuracy and reliability. Using this 
code, we have studied the matter and velocity density fields in 
the self-accelerating branch of the DGP model. We find good 
agreement of our matter power spectrum with that of previous 
work ||6ll l63ll . but the code enables us to go to even smaller 
scales and discover new interesting features. We also find that, 
like in the case of f(R) gravity Bill , the velocity field is more 
strongly affected by the modification to gravity law, and this 
provides a potentially powerful cosmological test of the the- 
ory of gravity. 

The structure of this paper is as follows. In § |II]we briefly 
describe the DGP model and the Vainshtein screening mech- 
anism, which lays down the essential theoretical background. 
In § |lll] we present the relevant field equations in the form 
as they appear in our numerical simulation code, and intro- 
duce the algorithm used to solve these equations as well as 
its implementation. In § [IV] we will show the various checks 
which we have performed to ensure that our code works accu- 
rately in different situations and aspects. We then run a suite 
of 36 simulations for the DGP model and its variants, and use 
these to study their predictions about the matter and velocity 
fields, especially on small scales which have been un-probed 
by previous studies: the results are given in § [V] Finally we 
summarise and outlook in § [VI] 

Throughout the paper we shall follow the metric convention 
(+, — , — , — ), and set c = 1 except in the expressions where c 
appears explicitly. 



II. THE DGP MODEL 



A. The model 



note that this branch of the DGP model suffers from the ghost 
instability (e.g. [680) and it cannot be considered as a physical 
model. Moreover, this model is already strongly disfavoured 
by observations, and therefore should be considered as a toy 
model to test our code. We will report results on more physical 
models as follow-up works in the future. 

The Friedman equation in the self-accelerating branch is 
given by 



H z = — 



H 8ttG 



-P, 



(1) 



in which r c is the cross-over length scale at which gravity be- 
comes five-dimensional, G is Newton's constant on the four- 
dimensional brane where matter fields live and p matter den- 
sity. To reproduce the observed cosmic acceleration one typi- 



cally needs r c 
ten as 



H . The Friedman equation can be rewrit- 



, , H(a) r— 
E(a) = — K -L = 



where a is the scale factor, a = 1 today, and 

1 „ 8nG 



PmO ■ 



(2) 



(3) 



The subscript indicates that the quantity is evaluated today. 

Under the quasi-static approximation, the Poisson equation 
and the equation for the scalar field are given by |j69| 



vV 



and 



3/3(a)a' 



■[(VV) 2 -(V i V j¥ »)(V i V^)] 



8nGa 2 
3/3{a) 



pS, 
(4) 



V 2 * = 4nGa 2 p6+-'V 2 tp, 



(5) 



where vf" is the gravitational potential, V is the spatial gradient 
operator and <5 is matter density contrast. The function (3(a) is 
defined b>Q 



/3 = 1- 2Hr, 



H 
3H 2 



which can be written as 



(6) 



(7) 



by 



^/firc (Qua 3 + fi rc ) 

If we linearise the equations, the Poission equation is given 



As an example of models that accommodate the Vainshtein 
mechanism, we consider the self-accelerating branch of the 
DGP braneworld model because this is the model where the 
Vainshtein mechanism has been studied in details (e.g., Ir65l - 
l67lo . This makes it possible to compare our new simulations 
with previous studies in the literature without the mesh refine- 
ments and demonstrate the power of the adaptive mesh tech- 
nique to study the nonlinear clustering. However, we should 



V 2 * = 4nGa 2 [l + —)p5. 



(8) 



Note that j3 is always negative so the growth of structure for- 
mation is suppressed in this model. 



We shall drop the time-dependence for j3 in the text hereafter for brevity. 
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B. Vainshtein mechanism 

The relevant field equation in a spherically symmetric con- 
figuration can be written in the following form 



(9) 



2r c 2 1 d 
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' 2 dtp~ 


3/3a 2 r 2 dr 




K dr) _ 


r 2 dr 


dr 



8ttG 

w 



6p m a 2 , 



in which 5p m = pS is the matter density perturbation and ip is 
the scalar field (i.e., the brane bending mode), and (3 is given 
in Eq. ((6]). Defining the mass enclosed in radius r as 



M(r) ee 4?r 
we can rewrite Eq. (O as 



Spmir 1 ] 



(10) 
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3/3 r 



dtp 
dr 



dr 



2 GM{r) 



^vM,(iD 



where for simplicity we have set a = 1, and is the Newto- 
nian acceleration caused by the mass M(r) at distance r from 
the centre. 

To further simplify the situation, let us assume that Sp m is 
constant within radius R and zero outside. Then Eq. (fTTT > has 
the physical solution 



dip 
dr 

for r > R and 

dip 
dr 
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(12) 



(13) 



for r < R. Here we have defined the Vainshtein radius r* by 

(14) 



„3 _ 8r c 2 r s 



* - 9/3 2 ' 

where r s = 2GM is the Schwarzschild radius. 

The fifth force is given by ^dtp/dr. So when r ^> r* we 
have 



1 dys 

2 d7 



1 

3/3 



(15) 



meaning that above the Vainshtein radius gravity is weakened 
(because f3 < for the self-accelerating branch of the DGP 
model). On the other hand, when r « r, we have 



(16) 



indicating that the fifth force gets suppressed well within the 
Vainshtein radius. 

As some useful examples, we have listed in the Table U the 
different radii of some objects. We could see that from objects 
as small as atoms to those as large as the Milky Way, the Vain- 
shtein radius is significantly larger than the physical size, and 
as a result the fifth force is negligible for those objects. 



TABLE I. The radii R (physical size), r s (Schwarzschild radius) and 
r» (Vainshtein radius) of some typical objects. Unit is meter. For es- 
timation we have used /3 = 2^2/3 and r c = H^ 1 ~ 4, OOOMpc ~ 
1.2 x 10 26 m. 



Object 


R 


r s 


r* 


Universe r 


- 1.2 x lO*" 


~ 4.5 x 10 25 r 


- 8.6 x 10 25 


Milky Way r 


- 0.9 x 10 21 


~ 2 x 10 15 


~ 3 x 10 22 


Sun 


- 0.7 x 10 9 


~ 3 x 10 3 ' 


- 3.5 x 10 18 


Earth 


~ 6 x 10 6 


~ 9 x 10" 3 


~ 5 x 10 16 


Atom 


- 5 X 10" 11 r 


- 1.8 x 10" 54 


~ 3 x 10" 1 



III. TV-BODY EQUATIONS AND ALGORITHM 

In this section we describe the A^-body equations in appro- 
priate units and their discretised versions which the ECOSMOG 
code solves to do cosmological simulations. 



A. TV-body equations 

The code unit used in our code is based on the superco- 
moving coordinates of 11331 17011 . which can be summarised as 
follows (tilded quantities are expressed in the code unit): 



x 



P = 



a 2 4> 



dt 



H 



pa 

dt 
7$' 



BH 

9 1 

c a ip 



(BH 



|2 ' 



(17) 



in which x is the physical coordinate, p c is the critical density 
today, Cl m the fractional energy density for matter today, v the 
particle velocity, <\> the gravitational potential and c the speed 
of light. In addition, B is the size of the simulation box in unit 
of /i _1 Mpc and Hq the Hubble expansion rate today in unit of 
lOO/i km/s/Mpc. Note that with these conventions the average 
matter density is p = 1. All these quantities are dimensionless. 

Using these quantities, the EOM of the brane-bending mode 
can be written as 



3/3a 4 



V 2 c£) - [ViVjip 







05-1*18) 



in which we have defined a new dimensionless 0(1) quantity 
H r c 1 



Rc. 



(19) 



For simplicity, from here on we shall ignore the tildes on all 
quantities in Eq. ( TPTI i. and unless otherwise stated these quan- 
tities are all in code unit where they appear. 

Eq. (fl~8l ) can be regarded as a second-order algebraic equa- 
tion for V 2 (^, which has two branches of solutions. To avoid 
undecidedness in which branch of solution (which may cause 
numerical problems in the code) to take, we first solve it ana- 
lytically and get 



1 



2(1 - w) 



-a ± ^a 2 +4(1 - w)E 



(20) 
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where we have defined 

_ 3/3a 4 



(21) 



E = (V 4 V^) 2 - w (VV) 2 + -B^ma (p - 1) , (22) 



in which we have followed [63] and used the operator-splitting 
trick to avoid the numerical problem due to imaginary square 
root, w is a constant coefficient which is chosen to be 1 / 3. 

Obviously, to get meaningful result in the limit one 
has to take the sign(a)-branch of Eq. (120b . where sign(a) is 
the signature of a or equivalently of (3, namely 



2(1 



+ T^yV '« 2 + 4 ( 1 - ™) S 



a 



(23) 



If a wrong branch is taken, then V 2 ^ ^ even when p — > 1, 
i.e., when the matter density field becomes uniform, and this 
is clearly a sign of inconsistency (indeed one of our code tests 
below is related to this observation). 

The equation for the standard Newtonian potential, on the 
other hand, is the same as in GR, and we have IT321 l33ll 



V 2 4> = -n m a(p- 1). 



(24) 



Once both <p and ip are solved, the total gravitational poten- 
tial in the DGP model, \I> = (f> + ip/2, follows, which can be 
differenced to find the total force under which particles move. 

I 



B. Discretisation 



Of course, computers could only handle a finite number of 
operations, and therefore Eqs. (T23ll24l i must be discretised be- 
fore being put into a code to be solved. This is the task for this 
subsection. 

Consider that we want to solve these equations on a three- 
dimensional mesh consisting of periodic cubic cells, with cell 
length h, and denote the value of the scalar field at the centre 
of the fc)-th cell by ipij,k (and similarly for other quanti- 
ties), then the discrete version of the field derivatives are 



^<P = ^ (fi+hj.k - <Pi-l,j,k) , 



(25) 



V 2 <<2 = ~2 (<Pi+l,j,k + fi-l,j,k - 2<Pi,j,k) , (26) 



V^VyV? = (<pi+l,j+l,k + <Pi-l,j-l,k — <Pi+l,j-l,k 



<Pi-l,j+l,k ) I 



(27) 



where we have assumed one dimension for simplicity for Wip 
and V 2 </?, the three-dimensional generalisations of which are 
trivial. 

Using these, it is straightforward to write down the discrete 
version of the Poisson equation as 



}?0 



'i+l,j,k + 4>i-l,j,k + <f>i,j+l,k + ^i.j-l.k + <f>i,j,k+l + <t>i,j,k-\ - §4>i,j,k) = ^m{Pi,j,k ~ 1) = -^m^i,j,k- (28) 



I 

Similarly, the EOM for the brane-bending mode can be writ- ten as an operator equation C h ipij t k = 0, with 
I 



£ Hi Pi,jM = j^(^Pi+l,j.k + <Pi-X,j,k + <Pi,j+l,k + <Pi,j-l,k + <Pi,j,k+l + Vi,j,k-1 — 6</9j,j,fe^ 



1 



2(1 - w) 



-a + o\/a 2 + 4(1 - w)T, itjtk 



(29) 



where the superscript h is used to label the level of the mesh (or equivalently the size of the cell of that level), and we have 

defined 
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1 — W 



h A 



(fi+l,j,k + fi-l,j,k — 2<Pi,j,k) + (ipi,j+l,k + <Pi,j-l,k — 2(P(.j,fc^ + (<Pi,j,k+l + <Pi,j,k-l — 2<Pi,j,k) 



"TjW (<fi+l,j,k + <Pi-l,j,k - 2<Pi,j,k) (<Pi,j+l,k + Vi,i-\,k - ^Vi,j,k) 

2 

"TjW (<fi+l,j,k + <Pi-l,j,k - 2<Pi,j,k) (<Pi,j,k+l + <Pij,k-l - 2<Pi,j,k) 

2 

"TjW (<Pi,j+l,k + <Pij-l,k - 2yi,j,fc) (<Pi,j,k+l + <Pi,j,k-l - 2yi,j,fc) 

1 / \2 

' g^4 yPi+l,3+hk + Vi— — l,fe — ¥>i+l,j-l,A — ^i-l,j + l,fc J 

+ ^4 (^i+lj,fe+l + Vi-l,j',fe-l - <Pi+l,j,k-l - <Pi-l,j,k+l) 

1 / \ 2 a 



(30) 



This equation can be solved using the Newton-Gauss-Seidel 
relaxation method, for which the code iterates to update the 
value of ¥>i,j,k in all cells, and at each iteration the field values 
changes as 



/i,ncw /i,old 

ri,j 3 k ^i,j,k 



h,old 



dip 



h, old 

6 a 

7? + 



4(1 -3w) 



(31) 



where we have 



/iV« 2 + 4(1 - wW^k 



Upi+l,j,k + <Pi-l,j,k + Pi,j+l,k + <Pi,j-l,k + <Pi,j,k+l + Pi,j,k-1 - 6</>i,j,fcJ • (32) 



Note that, in comparison with the corresponding equations 
in f(R) gravity 13211 . here the nonlinearity mainly exists in the 
term Sj^fc, which involves products of the derivatives. Such 
a special property of the Vainshtein screening makes it more 
difficult for the Newton-Gauss-Seidel relaxations to converge. 



C. Numerical implementation 

As mentioned above, RAMSES and therefore ECOSMOG are 
AMR codes. The code starts from a regular mesh which covers 
the whole simulation domain, which we call the domain grid, 
and adaptively refines the mesh if the effective number of dark 
matter particles in a given cell exceeds certain threshold Ar pc . 
In this way, it gives sufficient resolution in high matter den- 
sity regions and avoids wasting too much time on low density 
regions where resolution requirement is mild. The process of 
self-refinement goes on until finally all cells in the finest level 
do not satisfy the refinement criterion. Whenever a refinement 
is created, it is used to compute the force experienced by the 
particles which fall in its domain, and the boundary conditions 
can be set by using the values interpolated from coarser levels 



(see ll32"ll for more details). 

The essential features of the code that are relevant for mod- 
ified gravity simulations have been described in detail in ll32ll . 
Here we only present the most relevant part for our Vainshtein 
simulations. 

It is well known that, when using relaxation method to solve 
boundary value problems, the first few relaxations see the so- 
lution approach the true value quickly. The convergence then 
becomes slower as one gets closer to the true solution, making 
the method less efficient. This is because relaxation on the fine 
grid can only damp the short-wavelength Fourier components 
of the error, leaving behind the smooth long-wavelength com- 
ponents which cause a global poor convergence. The common 
method which is used to tackle this problem is the mutligrid 
technique lf73ll . The idea is simple: when the convergence rate 
becomes low on the fine grid, one coarsifi.es the discrete equa- 
tion by interpolation, moves it to the coarser level and solves it 
there. The long-wavelength components of the error then de- 
cays quickly on the coarse grids, giving more accurate solu- 
tion there, which can then be interpolated back to the fine grid. 
We follow the standard V-cycle in the arrangement of multi- 
grid: the field equation is always solved from the finest to the 



6 



coarsest grid, and then back. If convergence is not achieved 
(see below), further V-cycles are used. 

The multigrid method can be easily applied to computations 
based on a regular mesh, since the corresponding coarser grids 
are also regular and satisfy periodic boundary conditions. For 
refinements, which generally have irregular shapes, however, 
much effort needs to be devoted to designing the correspond- 
ing coarser meshes and setting appropriate boundary condi- 
tions for them. One of the important features of RAMSES and 
ECOSMOG is that it solves the discrete equation using multi- 
grid technique on both the domain grid and the refinements, 
which significantly speeds up the convergence. 

Consider the EOM for the brane-bending mode, which for 
simplicity can be written as 

C h y i ) = f h (33) 

on the fine level, where C is the nonlinear operator acting on 
tp h defined above. Note that although we have f h = for the 
DGP equation, we shall still keep f h in Eq. (1331 for the reason 
that will become clear soon. 

After a certain number of pre-smoothing Gauss-Seidel iter- 
ations on the fine grid, the convergence becomes slow because 
short-wavelength modes of the error have already decayed and 
long-wavelength modes are untouched. One arrives at an ap- 
proximate solution (p h on the fine level, which gives 

C h (<p h ) = f h ± f h . (34) 
Now consider what we call the residual equation, 

C h - C h (<p h ) = f h - f h = -d\ (35) 

where d h is the residual. After coarsifying and certain rear- 
rangements, we obtain the coarser-level equation as 

C H ( V H ) = C H (ll<p h ) - Ud h , (36) 

where 1Z is the restriction operator which compute the coarse- 
level values of quantities given their fine-level values by inter- 
polation. After a number of relaxation iterations on the coarser 
level, we find an approximate solution <p H , which has long- 
wavelength modes of the error further reduced, and therefore 
the fine-level solution can be corrected as 

where V is the prolongation operator which does the opposite 
of the restriction operator, again by interpolation. 
The truncation error r can be estimated as 

r h « C H (Tl0 h ) - TZC h (<p h ) , (38) 

and similarly for other levels. This could provide a stopping 
(convergence) criterion for the multigrid iterations, which in 
our case is 

\d h \<a\T h \, (39) 

in which a < 1/3 is some predefined constant, which is often 
set even smaller to be more conservative. Useful as it is, we 
find that this criterion is often not needed because it is easy to 
get \d h \ down to O (10 -10 ) within a few to a few tens of iter- 
ations, particularly on the refinements. This is typically much 
smaller than the truncation error. 



IV. CODE TESTS 

In this section we present the results of several tests of the 
ECOSMOG-V code, which are essential for us to be confident 
about its reliability. For this purpose, we have performed sev- 
eral tests as shown in Fig. Q] (see the figure caption for more 
technical details). 

A. Uniform density field 

The simplest test one can possibly do for the code is by con- 
sidering a homogeneous distribution of matter, in which case 
the scalar field ip remains a constant across the space, and its 
value can be arbitrary due to the shift symmetry in the EOM. 
If the code works correctly, then any inhomogeneity in <£>(x) 
should quickly disappear after a few Gauss-Seidel iterations. 

To check this, we set pi,j,k = 1 and use random values that 
follow a uniform distribution between —0.05 and 0.05 as the 
initial guess of (fi.j,k- These are shown in the upper left panel 
of Fig.[T]as the symbols, in which for clarity we have fixed the 
y and z coordinates. We then let the code do the Gauss-Seidel 
relaxation until the residual gets small enough, \d h \ < 10 -10 . 
After a few iterations this criterion is satisfied and the resulted 
ip(x) are described by the solid lines. We can see clearly that 
the final solution is constant in space. 

We have shown two different initial guesses as blue squares 
and red circles respectively, and tp(x) quickly converges to a 
constant, which is almost the same in both cases. 

B. One dimensional density field 

In one dimension we can check that (V 2 </?) 2 = (ViVjtp) 2 , 
which means that the nonlinear term in the DGP equation sim- 
ply vanishes. This provides another simple but very useful test 
of the code, because if there is something wrong in the nonlin- 
ear term in the code, then it is unlikely to vanish coincidentally 
and the solution for ip should be incorrect. 

Following ll32ll . we do two tests using one-dimensional den- 
sity field. The first test uses a sine-type field, specified by 

5{ x ) = --^—Ksm(27rx), (40) 

where if is a numerical coefficient which we take to be 

A' = 0.128tt 2 . (41) 

Because the nonlinear term in the EOM does not contribute in 
ID configurations, it can be easily checked that the analytical 
solution for this density field is 

ip(x) = 0.032 sin(27ra). (42) 

In the upper right panel of Fig. Q] we have presented this ana- 
lytical solution (blue solid curve) together with the numerical 
solution from the code (blue squares), which agree with each 
other very well. We have tried a few other numerical values of 
K and found same agreements. 
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FIG. 1. (Colour online) Code tests. Upper left panel: test using a uniform density field and random initial guess for </?(x) (§ HV Ab ; the symbols 
are the initial guesses and the solid lines are the final results of tp after relaxation. Upper right panel: tests using ID sine (blue squares) and 
Gaussian (red circles) density fields as described in Eqs. 1 1401143b respectively (§ 11 V Bb ; the solid curves are the analytical solutions. Lower 
panels: tests using a spherical overdensity with 8 = 23.77 (left) and S — 190.9 (right); the blue squares, red circles and black solid curves are 
respectively the code predictions on the domain grid, refinements and exact solutions obtained from numerical integration (§ IIV Cl and j HIVDb . 
In all tests we have assumed (3 = —5.0, a — 1.0; the simulation box size is 128/i _1 Mpc and the domain grid has 256 cells in each dimension. 



The second test uses a Gaussian-type density field, given by 



5{x) 



2Ja 



(a;-0.5) : 



x exp 



(x-0.5) : 



which corresponds to an exact analytic solution 
( \ rl"i ( (z-0.5) 2 



1 — a exp — 

y w~ 

Here J,a,w are constants which we take to be 

J = 0.02, a = 0.9999, w = 0.2, 



(43) 



(44) 



(45) 



as a numerical example. In the upper right panel of Fig. Q] we 
show this analytical solution (red solid curve) together with 



the numerical solution given by our code (red circles): again 
there is a very good match. 

Note that the code uses periodic boundary condition, which 
is exactly satisfied in the sine case. In the Gaussian case, the 
density field is not perfectly periodic, but S(x) is sufficiently 
suppressed by the exponential factor in Eq. ( 143b . and is close 
to zero at x = and x = 1, so that the periodic boundary 
condition approximately holds. 



C. Spherical overdensity 

In the above we have tested the code for ID density fields 
only. These tests show that the nonlinear term in the EOM do 
not create any problem, but they are not sufficient to prove that 
the code can solve this term correctly when it does contribute. 
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For the latter we need to do tests in three dimensions, of which 
the sphericall y sy mmetric cases are the simplest possibility. 

Following 16111 . here we assume that the matter density in- 
side the sphere is constant, i.e., we have a spherical overden- 
sity. In code unit, Eqs. (fl2l[T3T l can be written as 



dip 
dr 

for r < R and 

dip 
dr 



1L 



\ 1 8Cl m R 2 J _ 1 



9/3 2 



/ 1 | 8n m R2J R 3 i 



9/3 2 



(46) 



(47) 



for r > R, where r is the comoving coordinate scaled by the 
boxsize B, while R is the radius of the spherical overdensity 
also scaled by B. 5 is the overdensity. 

Given the value ip(r = 0), these equations can be integrated 
to find ip(r > 0) numerically. We do this for two different 
cases, respectively 8 = 23.77, R = 0.05 and S = 190.9, R = 
0.1, which correspond to the same mass in the sphere, and the 
results are shown as the black solid curves in the lower panels 
of Fig. Q] For comparison, solutions from our numerical code 
are shown as blue squares in those panels. 

We can see that in both cases the two solutions agree very 
well, especially on small r. Far from the centre, the agreement 
becomes less perfect since the numerical integration does not 
assume periodicity of the spherical overdensity, while the nu- 
merical code uses period boundary condition so that the spher- 
ical density sees its own images. Such an artificial discrep- 
ancy due to imposing periodic boundary conditions can also 
be seen in the point-mass test of l32ll . and it emphasises the 
importance of choosing large enough box sizes in cosmologi- 
cal simulations. 



D. Multilevel 

One of the most important features of our code compared 
with existing codes in the literature is that it enables adaptive 
mesh refinements to improve on the resolution. It is therefore 
crucial to check that the solver works correctly on the refine- 
ments. 

To check this, we have placed a refined region which covers 
the spherical overdensity in the above test. The density field 
on the refinement is chosen to the exactly the same as that on 
the domain grid, and we fix the boundary condition in the edge 
cells of the refinement by interpolating from the coarser cells 
in the domain grid. 

The results are shown as red circles in the lower panels of 
Fig. Q] where we can see that in the domain covered by the 
refinement they agree with both the solution from numerical 
integration (black solid curve) and the results on the domain 
grid (blue squares). We have also done tests using other values 
of S and found similar agreements. 

These tests show that the code actually works quite well on 
both the domain grid and refinements. In the next section, we 
will show results of cosmological simulations using this code. 



V. COSMOLOGICAL SIMULATIONS 

To study the cosmological behaviour of the DGP model and 
the Vainshtein mechanism, we have run a set of 36 TV-body 
simulations using the ECOSMOG-V code. In this section we 
show the results from these simulations. 



A. Simulation details 

As we have seen above, the self-accelerating branch of the 
DGP model shows a different background expansion history 
from that of ACDM, which itself can contribute to the dif- 
ferent matter clustering seen in these two models. Therefore, 
to see the effects of modified gravitational force and in par- 
ticular the Vainshtein mechanism more directly, it is better to 
compare the full DGP simulation with a corresponding New- 
tonian simulation with the DGP background: this so-called 
QCDM model ll6lll is not a rigorous theoretical model and we 
consider it only for comparisons. Also, we have simulated the 
linearised DGP model in which the nonlinear derivative cou- 
pling term is set to zero by hand so that the model has exactly 
the DGP background but particles feel an effective Newton's 
constant given by 



1 + 



(48) 



Because /3 < for the self-accelerating branch of the DGP 
model, the modified gravitational force is actually weaker than 
inGR. 

To make statistic averages, we run six realisations of sim- 
ulations for every model. The initial conditions are generated 
using the MPGRAFIC package at a = 0.02, from the same 
linear matter power spectrum but with different random seeds. 
Strictly speaking, even at Zi = 49.0 the matter clustering can 
be different in QCDM and DGP, but for the self-accelerating 
branch which we are studying here the difference in the linear 
matter power spectra is fairly small (~ 0.2%) and so we have 
neglected it. The fact that we use the same initial conditions 
for simulations of different models makes sure that the initial 
density fields have the same phases, and so any difference in 
the matter power spectra that we find at later times will be a 
direct consequence of the different dynamics and force laws 
between the models. 

The cosmological parameters we use in the simulations are 
taken from the best-fitting self-accelerating DGP model using 
the WMAP 5yr data lf74ll . To be more specific, we have 

{h, n s , In [lO 10 ^] , n b h 2 ,n c h 2 , n rc , a s } 
= {0.66, 0.998, 3.01, 0.0237, 0.0888, 0.138, 0.536}. 

Note that il rc is a derived parameter which makes the universe 
flat and erg is its current value for the linearised DGP model 
(because to compute it uses linear perturbation theory). Some 
of these parameters, such as f2 rc , f2 c , fi&, are used directly in 
the simulations, while all of them are used in generating the 
linear matter power spectra at z,. 
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FIG. 2. (Colour online) Visualisation of the density and velocity divergence fields from our high-resolution simulations with B = 100/i _1 Mpc 
and 256 3 particles on a domain grid which has 256 cells in each dimension. Upper panels: the density fields from one realisation of the QCDM 
(left), linearised DGP (middle) and full DGP (right) simulations, plotted on logarithmic scale; the bright (dark) regions have high (low) matter 
density. Lower panels: the velocity divergence fields from the same realisation of our QCDM (left), linearised DGP (middle) and full DGP 
(right) simulations; the colour bar indicates the values of the velocity divergence. All fields are plotted at a — 1.0. 



To better understand the effects of varying the force resolu- 
tion, we have simulated two different box sizes, respectively 
B = 100/i _1 Mpc and B = 200/i~ 1 Mpc. For both cases we 
use 256 3 dark matter particles and the domain grid has 256 
cells in each dimension. The cells are refined when the effec- 
tive number of particles inside them exceeds 9.0 (N pc = 9.0) 
and the final refinement level has 2 14 cells in each dimension 
if they were to cover the whole simulation box. The Gauss- 
Seidel relaxations stop when |e?' l | < e = 10~ 9 . 

Thanks to the efficient MPI parallelisation, the simulations 
are pretty fasj5 using 96 CPUs the low-resolution ones finish 



within 2 hours (wall-clock time) and the high-resolution ones 
complete in less than 6 hours. The difference is partly because 
the high-resolution simulations also have more refined time 
steps, and partly because of the stronger fluctuation of the den- 
sity field when constructed on a finer grid. Note that although 
the DGP simulations are in general more difficult to converge 
than the f(R) simulations (i.e., more tunings of the code are 
needed), if they do converge they do it much more quickly. 



2 Our DGP simulations take roughly 2-3 times CPU time as the correspond- 
ing QCDM runs. The simulations are expected to be slower if one increases 
the resolution, decreases the refinement criterion N pc or decreases the con- 
vergence criterion e. Also, in the self-accelerating branch of DGP model 
considered here, matter clusters less than in QCDM because the fifth force 



actually weakens gravity; in the normal-branch DGP, on the other hand, 
gravity is strengthened and matter cluster more strongly, which means there 
will be more refined grids and therefore the simulation takes longer to fin- 
ish. 
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B. Results 

1. Visualisation of density and velocity fields 

In Fig. [2] we have shown the density and velocity divergence 
fields in a slice of the simulation box. The velocity divergence 
field is defined as 

0(x) = lv • v(x), (49) 

in which H is the Hubble expansion rate and v is the veloc- 
ity field of the dark matter particles. To measure the velocity 
field and its divergence from the particle positions and veloc- 
ities we have used the Delaunay Tessellation Field Estima- 
tor (DTFE) code described in 07511 . The Delaunay tessellation 
has the advantage that the velocity divergence field computed 
in this way is volume-averaged, rather than mass-averaged. It 
also avoids the problem of empty cells including no particles, 
which can arise in direct assignment methods of measuring 
the velocity field IT76ll . 

Shown in the upper panels of Fig.|2]are the logarithmic den- 
sity field log p(x) for one of the 6 realisations of the QCDM 
(left), linearised DGP (middle) and full DGP (right) simula- 
tions. The colour scale is from p = 0.05 (black) to p = 100.0 
(white). The difference in the pattern and strength of cluster- 
ing is barely visible by eye, but one can still see the stronger 
clustering in the QCDM model (recall that the in DGP model 
gravity is weakened compared to QCDM), in particular in the 
cluster near the lower-right corner. A blink view of these pan- 
els also show that in QCDM the clusters are slight closer to 
each other, while the difference between linearised and full 
DGP is very small. 

Similar patterns can be seen in the lower panels of Fig. [2] in 
which we show the corresponding velocity divergence fields. 
The velocity divergence field is positive in low-density re- 
gions where matter flows from the central part and the flow 
becomes faster as it approaches clumps of matter. This trend is 
reversed near clusters and filaments, where the velocity diver- 
gence field becomes negative since matter flows inwards here. 
Inside clusters and filaments, the velocity divergence takes 
positive sign again, as noted by [ 76 J . Since gravity is strongest 
in the QCDM model, the matter flow is faster inside voids and 
near clusters, making the plot brighter in the low-density re- 
gions and the black regions near clusters thicker, though this is 
barely visible without blinking view (one can however see the 
difference in the clusters near the lower right corner of each 
panel). 

2. Matter and velocity divergence power spectra 

Given that the difference is so visibly small, it is more use- 
ful to look at statistical quantities, such as the power spectra, 
to quantify the matter clustering, 

Pss(k) = (|<5 k | 2 ), (50) 

with 

5 k = (2tt)~ 3/2 / <5(x) exp(-ik • x)d 3 x, (51) 



and 

Pee(k) = (|0 k | 2 >, (52) 

with 

9u = (27r) _3/2 J 0(x) exp(-ik • x)d 3 x, (53) 

for the velocity divergence field. In above (•) means ensemble 
average. 

In Fig. [3] we show the relative difference of Pgs (left panel) 
and Pqq of the full (blue symbols) and linearised (red sym- 
bols) simulations from that of the QCDM simulation. The left 
pan el agrees quite well with previous results such as those in 
16111 . but extends to smaller scales due to the AMR nature of 
our code. There are several notable features here: 

1. Both the linearised and the full DGP simulations agree 
with the linear perturbation theory prediction on length 
scales larger than k ~ O.lh Mpc -1 . This indicates that, 
at least for the WMAP5 best-fitting DGP model (the 
self-accelerating branch), the nonlinearity in Vainshtein 
mechanism becomes important at k ~ 0.1 h Mpc -1 , 
exactly where the nonlinearity in the underlying density 
field starts to make linear perturbation theory invalid. 

2. On scales smaller than k ~ 0.1/iMpc -1 , the Vainshtein 
effect can strongly suppress the deviation from QCDM 
compared to the linearised simulation (where it is ab- 
sent). This agrees with the result of ll6lll . 

3. The AMR property of our code enables us to measure 
the matter power spectra on significantly smaller scales 
than that of ll6lll . Indeed, the simulations show that the 
deviation from QCDM in P$s decays to zero on small 
scales, in agreement with the halo model prediction of 
& 

4. The low-resolution simulations, with a larger box size, 
show better agreement with linear perturbation theory 
on large scales as expected. Box size is an equally, if not 
more, important issue in our simulations, and should not 
be chosen to be too small as a compromise to achieve 
better resolutions. This reflects the importance of AMR 
from a different angle. 

In the framework of GR, fitting formulae such as the Halofit 
have been developed to provide a mapping between the linear 
and nonlinear matter power spectra nfm . But due to its sim- 
plicity, it is often misused to derive the nonlinear power spec- 
trum in modified gravity models IT79ll . Fig.l3lclearly shows that 
the Halofit overestimates the deviation of DGP from QCDM 
on small scales. This is hardly surprising because Halofit does 
not incorporate the effect of the Vainshtein mechanism which 
helps to recover GR on small scales. The Vainshtein mecha- 
nism is effective once the density perturbations become non- 
linear l69ll . and brings the matter power spectrum back to the 
QCDM prediction on small scales. Based on this observation, 
118011 has suggested a simple way to modify the Halofit to fit the 
nonlinear power spectrum in modified gravity models, which 
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FIG. 3. (Colour online) The relative differences of the matter (left panel) and velocity divergence (right panel) power spectra of the full (blue 
symbols) and linearised (red symbols) DGP simulation from that of the QCDM simulation, at a — 1.0. The results are obtained by averaging 
6 realisations for the high- and low-resolution simulations respectively (see the texts in the figure for more details). The velocity divergence is 
measured from a Delaunay tessellation of the particle distribution. We also show the predictions of linear theory (black dashed curves), Halofit 
(black dotted curve) and PPF (black solid curve). 



-cos rr 




k (h Mpc ') k (h Mpc ') 

FIG. 4. (Colour online) Same as Fig. |3]but for three different output times - a = 0.3 (dotted curves), a — 0.5 (dashed curves) and a = 1.0 
(solid curves) - to see the time evolution of the pattern. Blue squares are from the full DGP simulations and red circles from the linearised 
simulations. For clarify we have only shown the results from the B = 100/i~ Mpc simulations and error bars are not plotted. The smooth 
black curves in the left (right) panel are predictions of PPF (linear perturbation theory). 
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is called the Parameterised-Post Friedman (PPF) fit. The PPF 
matter power spectrum is given by 

„ ,, \ Ps8,non-GR .QCDM 

ns{ *> Z) ~ l + c nl (z)£ 2 (M) 

(54) 

where z) measures the nonlinearity of the density pertur- 
bations 

u3 

X 2 {k,z) = —P S 5,lin, (55) 

Pss,\in is the linear matter power spectrum for the DGP model 
and Paa.non-GR is the power spectrum in the linearised model 
without the Vainshtein mechanism, which can be obtained by 
simply applying Halofit to the linear power spectrum for the 
DGP model. Pss, qcdm is the power spectrum in the QCDM 
model. As the nonlinearity becomes large S > 1, this for- 
mula ensures that the power spectrum in the full DGP model 
approaches the one in the QCDM model. In Ref. 18111 . the free 
parameter c n \(z) in the PPF fit was estimated using the pertur- 
bation theory as c„i(0) = 0.3, c n \(l) = 0.34 and c n i(2.33) = 
0.36. Using these predictions for c n \(z), we found that the PPF 
fit recovers the iV-body simulation results very well (Figs. [3] 
and0J. This demonstrates that the Vainshtein mechanism is 
working as expected in our full DGP simulations. 

The right panel of Fig.[3]is the same as the left panel, but for 
the velocity divergence power spectra. Here, the behaviour is 
very different. In particular, we see that the Vainshtein mecha- 
nism already has effects on sales as large as k ~ 0.05/iMpc - , 
similar to what we have seen in our f(R) gravity simulations 
ll3"Tll . Furthermore, the deviation from QCDM on large scales 
(~ 30 — 35%) is much larger than what we see in P$s (~ 13% 
as is shown in the left panel). This means that observationally 
the velocity field, though more challenging to measure, can be 
more interesting in the tests of gravity. 

3. Time evolution of the power spectra 

Finally, Fig. |4] shows the time evolutions of the matter (left 
panel) and velocity divergence (right panel) power spectra. In 
particular, we show AP SS / Pss,qcbm and AP^/P^^qcdm 
at three different epochs, with a = 1.0 (solid curves), a = 0.5 
(dashed curves) and a = 0.3 (dotted curves). Results for the 
full (linearised) simulations are shown in blue squares (red 
circles) as before. 

The left panel of Fig. [4] shows that deviations from QCDM 
are smaller at earlier times, which is because \j3\ is larger and 
G e g is closer to Gn at higher redshifts. The Vainshtein effect 
is also weaker at the early times because the nonlinear term in 
the DGP equation is proportional to \f3\~ 1 and therefore sup- 
pressed. However, the general observations that we have for 
a = 1.0, notably that the Vainshtein mechanism suppresses 
the deviations from QCDM on smaller scales and eventually 
brings things back to GR, still apply here. The evolution pat- 
tern can be understood as follows: 

1 . Note that the deviation from QCDM essentially freezes 
on scales smaller than k ~ 0.7/iMpc -1 after a = 0.3, 



which shows that those scales become below the Vain- 
shtein radius so that they feel the standard gravity. 

2. The very large scales, on the other hand, are beyond the 
Vainshtein radius until today and their evolution can be 
approximately described by linear perturbation theory. 

3. Going from very large scales to intermediate scales, the 
nonlinearity of the matter density field first kicks in and 
makes the full DGP simulation behave similarly to the 
linearised simulation (i.e., AP$s /Pss initially decreases 
as k increases). Then, the Vainshtein mechanism takes 
over and brings things back to the small-scale behaviour 
described above. This explains why APgs/Pss appears 
to have a valley whose position depends on the cosmic 
epoch and shifts to larger scales with time. 

The velocity divergence power spectrum, on the other hand, 
shows very different properties from the matter power spec- 
trum: 

1. Even on length scales as large as k ~ 0.1/iMpc -1 , Pgg 
in the full DGP simulations does not agree with linear- 
theory prediction, though in the linearised DGP simu- 
lations it does agree with linear theory on those scales. 
This implies that the Vainshtein mechanism affects the 
velocity divergence even on such large scales, and it is 
in contrast to the matter power spectrum, which agrees 
with linear theory predictions in both the full and the 
linearised DGP simulations on those scales. This has an 
important implication in the measurement of the growth 
rate using redshift distortions. We need to carefully take 
into account the Vainshtein effect even on large scales 
to extract the linear growth rate accurately. 

2. On small scales, the Vainshtein mechanism suppresses 
the deviation from QCDM,a but the velocity divergence 
power spectrum does not approach to the QCDM re- 
sult even at k ~ 10/iMpc -1 whereas the matter power 
spectrum does. This implies that velocities hold a key to 
distinguish modified gravity models with the Vainshtein 
mechanism from GR on small scales. 

We can also understand the pattern for APgg/Pgg. Taking 
the linearised DGP simulation as an example, for which New- 
ton's constant is scaled by 1 + 1/(3/3) on all scales, so that in 
linear theory one would expect APgg / Pqq to be a constant at 
least on large scales, as we could see at early times (a = 0.3). 
On small scales, the nonlinearity transfers power between dif- 
ferent scales, causing the deviation from scale-independence. 
When clusters start to form, the 8 field changes sign from the 
outside to the inside of halos, and its magnitude can become 
smaller. As gravity is stronger in the QCDM model, at a given 
time it has stronger clustering and more compact clusters, and 
so the scale at which 8 changes sign is also smaller - this leads 
to, on average, an increase in \APgg\/ Pgg on the typical halo- 
formation scale and a decrease on slightly larger scale, which 
gives rise to the peak we see in the plot of APgg/Pgg. The 
peak shifts to larger scales because the halo-formation scale 
increases in time. 
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VI. SUMMARY AND OUTLOOK 

To summarise, in this paper we have reported the develop- 
ment and initial results of a new A^-body code which is suit- 
able for cosmological simulations in modified gravity theories 
which restore GR in high-density regions by nonlinear deriva- 
tive couplings, such as the DGP and Galileon gravity mod- 
els. This code is based on our ECOSMOG code [3211 . which 
was originally developed for f(R) gravity simulations and 
later generalised to simulate generic modified gravity theories 
which restore GR using nonlinear self-interaction potentials 
1361 l37ll . However, the solver for the nonlinear partial differ- 
ential equation is different in both its algorithm and its im- 
plementations, making it essentially a new and independent 
code. 

The code, which we call ECOSMOG-V with the V standing 
for 'Vainshtein', has a few distinct features compared with ex- 
isting codes in the literature which tackle the same problem: 

1 . It solves the scalar field on a mesh using nonlinear re- 
laxation, which has good convergence properties ll63l . 
The density field on the mesh is obtained by assigning 
particles following the triangular-shaped-cloud scheme, 
and there is no need to pre-smooth it to achieve conver- 
gence. 

2. The mesh can be adaptively refined in high-density re- 
gions to achieve higher resolution and accuracy there, 
without affecting the overall performance. There is no 
need to smooth the density field on the refinements ei- 
ther. 

3. It is efficiently parallelised using MPI, which makes the 
simulations fast. This is important for modified gravity 
simulations which generally take much longer than cor- 
responding GR simulations, and makes this code suit- 
able for large and high-resolution cosmological simula- 
tions which are essential for us to fully understand the 
small-scale behaviour of such models, and to better ex- 
plore the future high-precision observational data on the 
large-scale structures. 

We have made various tests of the code to make sure that it 
gives accurate solutions to simple situations such as uniform 
density field, ID (sine and Gaussian) density fields and spher- 
ical overdensites. Because the most important feature of this 
code is the AMR, we also tested the code on refinements and 
found that it works correctly. Furthermore, our cosmologic al 
simulations predict the same behaviour of the matter power 
spectrum as that discovered in previous work I16lll63l . 



We have used our code to run a set of 36 simulations for the 
QCDM, linearised DGP and full DGP models. The AMR na- 
ture enables us to probe scales smaller than those reached by 
previous simulations, and confirmed the halo-model predic- 
tion of 1 77] that the full DGP matter power spectrum should 
go back to the QCDM result on small enough scales. 

We have also studied the velocity field in DGP gravity, and 
found that it is more strongly affected by the nonlinearity in 
both the underlying matter field and the Vainshtein mecha- 
nism. Even on large scales where the matter power spectra 
for the full and linearised DGP simulations agree quite well, 
we find noticeable difference between the velocity divergence 
power spectra of the two. This trend starts as early as z = 1, 
and by z = the deviation of Pee from the corresponding 
QCDM result can be as large as more than 30%, compared 
with the ~ 13% deviation in P$s- This suggests that large- 
scale velocity fields can be a good probe of modified gravity 
theories such as the DGP. 

Of course, the DGP model is disfavoured by various ob- 
servational tests, but the beauty in the idea of Vainshtein 
mechanism has motivated more general models which recover 
GR using nonlinear derivative couplings, notably the Galileon 
models. Recent studies of f58il60ll have largely advanced our 
knowledge in the linear perturbation behaviour of these mod- 
els and found reasonable fits to the CMB and background ex- 
pansion data, but still leave the open question about the re- 
gion of validity of linear theory. To fully answer this question 
is crucial for cosmological tests of the Galileon model, espe- 
cially in light of the coming high-precision observational data, 
but this can only be achieved by studying the nonlinear regime 
of the structure formation numerically. Our code is useful in 
this context, and could hopefully bring us one important step 
forward. 
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